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The Kondo lattice model describes a quantum phase transition between the antiferromagnetic 
state and heavy-fermion states. Applying the dual-fermion approach, we explore possible super¬ 
conductivities emerging due to the critical antiferromagnetic fluctuations. The d-wave pairing is 
found to be the leading instability only in the weak-coupling regime. As the coupling is increased, 
we observe a change of the pairing symmetry into a p-wave spin-singlet pairing. The competing 
superconductivities are ascribed to crossover between small and large Fermi surfaces, which occurs 
with formation of heavy quasiparticles. 


Heavy-fermion superconductivities in rare-earth and 
actinide materials commonly appear near antiferromag¬ 
netic (AFM) or ferromagnetic quantum critical point 
(QCP). This situation resembles high-Tc cuprates, where 
d-wave superconductivity (d-SC) with the gap symmetry 
(j)k oc cos kx — cos ky is realized due to AFM fluctuations 
of q = (7r,7r) = Q. This similarity suggests a possibil¬ 
ity of the d-SC also in heavy-fermion systems near AFM 
QCP. However, since the microscopic origin of the mag¬ 
netism is different between d- and /-electron systems, 
pairing interactions derived from the magnetic fluctua¬ 
tions may be dissimilar, resulting in a different type of 
superconductivity in /-electron systems. 

Motivated by the above idea, we consider the follow¬ 
ing simplified situation to focus on a difference in the mi¬ 
croscopic interactions. We first suppose a tight-binding 
band on the square lattice, = —2{coskx + cosky), 
which exhibits a perfect nesting at gr = Q at half-filling. 
In the case of the cuprates, the Hubbard interaction U 
gives rise to the AFM long-range order, and carrier dop¬ 
ing leads to the d-SC. In heavy-fermion systems, on the 
other hand, the nesting of the conduction band leads to 
the AFM ordering of localized /-electrons through the 
RKKY interaction. Our question now is whether the d- 
SC is realized in this case as well by tuning the system 
toward QCP, and if not, which type of pairing is favored? 

For this purpose, we explore possible unconventional 
superconductivities in the two-dimensional Kondo lattice 
model. The Hamiltonian reads 

^ ^ J ^ ^ ( 1 ) 

ka i 

where Si = ^he number of lat¬ 

tice sites, N, is taken to be Y = with L = 32,64,128. 
Various numerical methods have been applied to the 
Kondo lattice with central attention on the magnetic 
properties [1-6]. Further elaborate calculations are nec¬ 
essary to elucidate pairing fluctuations lying near the 
magnetic instability [7]. 

In order to address superconductivity in the Kondo 
lattice, we need an approximate method which, at least. 


captures Kondo physics and AFM fluctuations around 
QCP. The dynamical mean-field theory (DMFT) takes 
full account of local correlations [8] , and hence of Kondo 
physics. Indeed, formation of heavy quasiparticles was 
traced at finite temperatures [9], and a quantum phase 
transition between the Kondo paramagnet and the AFM 
state was derived [10, 11]. We work on an extension of the 
DMFT to describe influence of critical AFM fluctuations 
on heavy quasiparticles. 

The dual-fermion approach provides a way to per¬ 
form diagrammatic expansion around the DMFT [12, 13]. 
In particular, inclusion of ladder-type diagrams as in 
the fluctuation exchange (FLEX) approximation describe 
long-range correlations, which lead to paramagnon exci¬ 
tations [14, 15], the d-SC in the Hubbard model [16], 
and the correct critical exponents in the Falikov-Kimball 
model [17]. Further advantage of the dual-fermion 
method as compared with related theories [18-21] is its 
applicability: The framework is easily applied to the 
Kondo lattice model as presented below. 

Formalism. In the path integral representation us¬ 
ing Grassmann numbers, the partition function of the 
model (1) is written as Z = f Yli[F{c*Ci)VSi]e~'^ with 
the action 

S = ~ ^uj)ctjka^ujka- ( 2 ) 

i Ldka 

Here is the local part defined by 

^imp ^ +11- Au)c+^Cu;ia + JSi ■ Si- (3) 

UJCT 

The dual-fermion transformation is performed on the 
kinetic-energy term described by the second term of S. 
On the other hand, the difference between the Kondo lat¬ 
tice model and the Hubbard model lies in the local part 
^imp xherefore, the procedure deriving the dual-lattice 
action in the Hubbard model can be directly applied to 
the Kondo lattice, leaving the difference in the local in¬ 
teractions. Introducing dual fermions d and formally in¬ 
tegrating out the variables c and S at each site, we arrive 
at the lattice model Z = /T V{dldi)e ^ which is writ- 
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FIG. 1. Diagrams for (a) the dual self-energy E and (b) the 
pairing vertex F^^. The lines and squares represent G and 7 , 
respectively. 
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Information of the local interaction enters in and 
7 . The propagator & is defined by = {g~^ + 

- Ck)-^ - with g^ = -(c,^iaC*^a)imp being the 
local Green’s function in the system described by 5]™^. 
The interaction coefficient 7 corresponds to the vertex 
part in the impurity system: 71234 = [(ciC 2 C 3 c 2 )inip — 
^i 5 ' 2 (^ 14^^23 -^^ 13 ^^ 24 )]/(T^i^ 25 ' 35 ' 4 ) where we used abbre¬ 
viations such as 1 = (co’i,cri). In deriving Eq. (4), we 
neglected many-body interactions which involve more 
than three particles. Validity of this approximation 
was confirmed numerically [14] and also supported in 
terms of 1/d expansion [16]. We compute g^j and 7 us¬ 
ing the continuous-time quantum Monte Carlo (QMC) 
method [ 22 ] applied to the Kondo impurity model [23]. 

Our task is now to evaluate the dual self-energy 
which is connected to the original Green’s function Gujk 
by the exact formula [12, 13], Gujk = [{guj^guj'^ujkguj)~^ + 
Auj — €k]~^• We take FLEX-type diagrams in the particle- 
hole channel expressed in Eig. 1(a) to take account of crit¬ 
ical AEM fluctuations [14-16] [24]. After J^^k is obtained, 
we update the hybridization function and solve again 
the effective impurity problem until the self-consistency 
condition G^jk = 0 is fulfilled. 

QCP. We first identify the quantum phase transition 
between Kondo insulator and AEM insulator at half¬ 
filling. According to a lattice QMC calculation [ 1 ], the 
critical coupling is estimated at = 1.45 ± 0.05. 

Comparing with this value, we shall check accuracy of 
our approximation. 

In the dual-fermion approach, the AEM fluctuations 
may be observed through eigenvalues of the matrix 
Aujuj' = -(l/V)y;fcGa;fcGa;.fc+Q7X' [16]- Here 
denotes the spin channel of the local vertex, 7 ®^, = 
The transition takes place when the 
largest eigenvalue Aafm exceeds 1 . Because of the 
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FIG. 2. (Golor online) (a) The leading eigenvalue Aafm for 
the AFM fluctuations at half-filling, n = 1. The dashed 
lines indicate the critical behavior 1 — Aafm oc (b) 

An intensity plot of Aafm in the J-T plane. The points indi¬ 
cate QGP estimated by QMG, DMFT, and our approximation 
(DF). The dashed line shows the phase boundary computed 
in DMFT. 


Mermin-Wagner theorem [25], the AFM transition is for¬ 
bidden at T 7^ 0 and Aafm follows the critical behavior 
1 — Aafm oc when the ground state is AFM [26, 27]. 
On the other hand, if the AFM is suppressed by the 
Kondo effect, Aafm becomes constant at low tempera¬ 
tures. To distinguish these behaviors, we plot 1 — Aafm 
as a function of l/T in Fig. 2(a). These plots demon¬ 
strate that the ground state is AFM for J < 1.2 and 
paramagnetic for J > 1.4. The data at J = 1.3 show 
a clear size dependence, which is common to the AFM 
ground state. Therefore, we conclude that J = 1.3 is 
in the AFM region, and the quantum phase transition is 
estimated at = 1.35 ± 0.05. This value agrees with 
the QMC result within 10 %, and shows considerable im¬ 
provement from the DMFT result = 2.18. These 

estimations of Jc are summarized in Fig. 2(b) together 
with an intensity plot of the AFM fluctuations. The crit¬ 
ical region, e.g., 1 — Aafm ^ 10 “^, forms a dome which 
is similar to the phase boundary in the DMFT. 

Superconductivity. We now compute superconduct¬ 
ing instabilities near half-filling. We solve the lin¬ 
earized Bethe-Salpeter equation for pairing correlations, 
-(T/N) Yk' GuG-uTll,<i>k' = Asc</>fe, where k = 7,fc)- 
For the irreducible vertex part F, we take account of 
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strong-intensity regions with opposite signs. It indicates 
that those pairings are due to AFM fluctuations. The 
angular dependence for Big and are expressed as 

oc co^kx — cos ky^ (5) 

(/)p(xsin% {^ = x,y), (6) 

which may be referred to as d-wave and p-wave, respec¬ 
tively. The transition temperature determined from the 
condition Asc = 1 is shown in Fig. 4(b). The d-SC occurs 
in the weak-coupling regime J < 0.9, while the p-wave 
superconductivity has a higher transition temperature in 
J > 0.9. We observed a similar behavior in the wide- 
doping range of n > 0.80. 

It is noteworthy here that the competing d-wave and 
p-wave superconductivities have been discussed in a con¬ 
text of CeRhIns and CeCu 2 Si 2 [29]. By a phenomeno¬ 
logical treatment of AFM QCP, properties of the p-wave 
singlet superconductivity was investigated. 

As J is increased further beyond Jc, on the other hand, 
we observed the leading 5-wave fluctuations. Figure 3(b) 
shows Asc and (j)ujQk for J = 1.8 and n = 0.92. The 
eigenfunction of the 5-wave fluctuations has the momen¬ 
tum dependence expressed by 


FIG. 3. (Color online) Temperature dependence of the leading 
eigenvalues Asc for pairing fluctuations. The intensity plots 
show the corresponding eigenfunctions (j)ujQk in the momentum 
space for the three largest fluctuations. The arrows indicate 
the nesting vector Q = (7r,7r). (a) J = 0.8, n — 0.84 and 
(b) J = 1.8, n = 0.92. The system size is A = 32^. 

ladder-type vertex corrections in Fig. 1(b) to describe 
pairing interactions mediated by AFM fluctuations [16]. 
As in the case of AFM ordering, the superconducting 
transition can be detected by the condition Asc = 1 - 
We computed the leading eigenvalues Asc for all types of 
pairing symmetry, i.e., 2 spin channels (singlet/triplet) 
and 5 spatial symmetries classified according to the irre¬ 
ducible representations of the point group D 4 . 

Figure 3(a) shows temperature dependence of Asc for 
n = 0.84 in the weak-coupling regime, J = 0.8 < Jc. 
The leading instability is the spin-singlet pairing with Big 
symmetry as expected. However, other symmetries in 
the spin-singlet channel show comparable enhancement 
in contrast to the Hubbard model [16]. In particular, the 
Eu pairing is nearly degenerate with Big and becomes the 
leading instability in J = 1.0 as also plotted in Fig. 3(a). 
We note that the Eu pairing is an odd function in the 
frequency domain as well as in the momentum space [28] . 
The momentum dependence of the eigenfunctions (j)ujok 
at the lowest Matsubara frequency ujo = ttT is plotted in 
Fig. 3(a) for the three leading fluctuations. A common 
feature is found: The nesting vector Q = (tt, tt) connects 


oc cos kx + cos ky^ (7) 

which may be referred to as an extended-<s wave. The 
enhanced 5-wave pairing was observed in the low-doped 
and strong-coupling regime, n > 0.88 and J ^ Jc- How¬ 
ever, we did not find a transition to the 5 -wave super¬ 
conductivity, i.e., a parameter set which gives Asc = 1, 
or could not reach low enough temperatures because of 
a bad convergence of the bath function A^;. 

In order to identify the driving force of the 5 -wave fluc¬ 
tuations, we performed DMFT calculations without dual 
fermions as well. In this case, no enhancement of the 
fluctuations was observed within temperatures we have 
reached [30]. Therefore, we conclude that spatial corre¬ 
lations are relevant for enhancement of the 5-wave fluc¬ 
tuations. An 5 -wave superconductivity was recently re¬ 
ported in the periodic Anderson model [31], and its con¬ 
nection with our results is left for future investigations. 

Quasiparticle band. We discuss quasiparticle proper¬ 
ties to reveal the origin of the competing superconduct¬ 
ing fluctuations. Two contrasting spectra A{k^uj) = 
—^:~Am.G^^k are shown in Fig. 4(a). Here, we used 
the Fade approximation for analytic continuations. The 
non-interacting dispersion Ck can be traced in J = 0.8, 
while the spectrum exhibits a distinct hybridization gap 
in J = 2.6 as a consequence of the Kondo screening. The 
formation of quasiparticles can be quantified in terms of 
the topology of the Fermi surface (Fermi line in two di¬ 
mensions) [9]. To this end, we evaluated the area nps 




























4 




(71,0) (0,0) (7l,7l) (71,0) (71,0) (0,0) (7l,7l) (71,0) 

(71/2,71/2) (71/2,71/2) 


77 = 0.84 TZps 



0 0.5 1 1.5 2 2.5 3 

J 

FIG. 4. (Color online) (a) The single-particle excitation spec¬ 
trum A{k,uj) for n = 0.84 and T = 0.056. (b) Supercon¬ 
ducting transition temperatures of singlet-Big (circle) and 
singlet-(square) symmetries. The intensity plot shows the 
Fermi-surface area nrs defined in Eq. (8). The Fermi-surface 
structure is plotted for three representative parameters. 

surrounded by the Fermi line: 

2 

'^FS ~ ^ ^ R'6^u;=0,fc)5 (8) 

k 

where 0 is the step function. A crossover between the 
non-interacting band (nps = '^) and the heavy-fermion 
band (nps = n + 1) is thus visualized in the J-T plane 
in Fig. 4(b). It is worth notice that the crossover region 
coincides with the coupling region where the change of 
the pairing symmetry is observed. 

In what follows, we discuss the origin of the compet¬ 
ing superconducting fluctuations in connection with the 
crossover of the Fermi-surface structure. As we have 
pointed out, both d-wave and p-wave pairings are in¬ 
duced by AFM fluctuations. The d-wave is in partic¬ 
ular favorable in the weak-coupling regime because the 
intensities of (j)ujQk lie around the van Hove points, i.e. 
ki = (tt, 0 ). The advantage of the d-SC, however, dimin¬ 
ishes as the hybridization band is formed. The p-wave 
superconductivity instead emerges in the crossover re¬ 
gion between small and large Fermi surfaces as shown 
in Fig. 4(b). In this region, low-energy excitations exist 
around k 2 = ( 7 r/ 2 , 7 r/ 2 ) rather than around the van Hove 
points ki. Therefore, the scattering between k 2 and —k 2 
gives a major contribution to pairing interactions to favor 
the p-wave symmetry. Furthermore, dynamical nature of 


quasiparticle interactions is relevant for a realization of 
odd-frequency pairings [28, 29]. In the present formalism, 
the vertex 7 plays its role, and it indeed has a strong fre¬ 
quency dependence in the heavy-fermion regime. 

In summary, the d-SC is not particularly favored in 
the Kondo lattice even though the nesting at gr = (tt, tt) 
gives rise to critical AFM fluctuations. The key fea¬ 
ture which yields differences with the Hubbard model 
is the Fermi-surface crossover: The nesting of the non¬ 
interacting band governs the magnetic instability, while 
the hybridization band is formed over the evolution of 
heavy quasiparticles. The “incomplete” quasiparticles in 
the crossover regime undergo the critical AFM fluctua¬ 
tions, resulting in emergent quantum states such as the 
p-wave spin-singlet superconductivity. 

The author thanks H. Kusunose for encouraging dis¬ 
cussions. This work was supported by JSPS KAKENHI 
Grant No. 26800172. A part of the computations was 
performed in the ISSP Supercomputer Center, the Uni¬ 
versity of Tokyo. 


[1] E. E. Assaad, Phys. Rev. Lett. 83, 796 (1999). 

[2] H. Watanabe and M. Ogata, Phys. Rev. Lett. 99 , 136401 
(2007). 

[3] L. G. Martin and E. E. Assaad, Phys. Rev. Lett. 101 , 
066404 (2008). 

[4] T. Misawa, J. Yoshitake, and Y. Motome, Phys. Rev. 
Lett. 110 , 246401 (2013). 

[5] G. Li, arXiv:1309.0156 (2013). 

[ 6 ] R. Peters and N. Kawakami, arXiv: 1504.07338 (2015). 

[7] M. Z. Asadzadeh, E. Becca, and M. Eabrizio, Phys. Rev. 
B 87 , 205144 (2013). 

[ 8 ] A. Georges, G. Kotliar, W. Kranth, and M. J. Rozen- 
berg. Rev. Mod. Phys. 68 , 13 (1996). 

[9] J. Otsuki, H. Kusunose, and Y. Kuramoto, Phys. Rev. 
Lett. 102, 017202 (2009). 

[10] R. Peters and T. Pruschke, Phys. Rev. B 76 , 245101 
(2007). 

[11] J. Otsuki, H. Kusunose, and Y. Kuramoto, J. Phys. Soc. 
Jpn. 78 , 034719 (2009). 

[12] A. N. Rubtsov, M. 1. Katsnelson, and A. 1. Lichtenstein, 
Phys. Rev. B 77 , 033101 (2008). 

[13] A. N. Rubtsov, M. 1. Katsnelson, A. 1. Lichtenstein, and 
A. Georges, Phys. Rev. B 79 , 045133 (2009). 

[14] H. Hafermann, G. Li, A. N. Rubtsov, M. 1. Katsnelson, 
A. 1. Lichtenstein, and H. Monien, Phys. Rev. Lett. 102 , 
206401 (2009). 

[15] H. Hafermann, Numerical Approaches to Spatial Correla¬ 
tions in Strongly Interacting Fermion Systems (Guvillier 
Verlag, Gottingen, 2010). 

[16] J. Otsuki, H. Hafermann, and A. 1. Lichtenstein, Phys. 
Rev. B 90 , 235132 (2014). 

[17] A. E. Antipov, E. Gull, and S. Kirchner, Phys. Rev. 
Lett. 112 , 226401 (2014). 

[18] H. Kusunose, J. Phys. Soc. Jpn. 75, 054713 (2006). 

[19] A. Toschi, A. A. Katanin, and K. Held, Phys. Rev. B 
75 , 045118 (2007). 

[20] G. Slezak, M. Jarrell, T. Maier, and J. Deisz, J. Phys.: 













5 


Condens. Matter 21, 435604 (2009). 

[ 21 ] C. Taranto, S. Andergassen, J. Bauer, K. Held, 
A. Katanin, W. Metzner, G. Rohringer, and A. Toschi, 
Phys. Rev. Lett. 112 , 196402 (2014). 

[22] E. Gull, A. J. Minis, A. I. Lichtenstein, A. N. Rubtsov, 
M. Troyer, and P. Werner, Rev. Mod. Phys. 83 , 349 
( 2011 ). 

[23] J. Otsukl, H. Kusunose, P. Werner, and Y. Kuramoto, 
J. Phys. Soc. Jpn. 76 , 114707 (2007). 

[24] The ladder approximation Is sufficient for our purpose 
of revealing pairing symmetries Induced by the critical 
AFM fluctuations. More complicated diagrams need to 
be Included to address, e.g., feedback from pairing fluc¬ 
tuations, which Is beyond the scope of this paper. 

[25] N. D. Mermln and H. Wagner, Phys. Rev. Lett. 17 , 1133 
(1966). 


[26] S. Ghakravarty, B. 1. Halperln, and D. R. Nelson, Phys. 
Rev. Lett. 60 , 1057 (1988). 

[27] P. Hasenfratz and F. Nledermayer, Physics Letters B 
268 , 231 (1991). 

[28] A. Balatsky and E. Abrahams, Phys. Rev. B 45 , 13125 
(1992). 

[29] Y. Fuseya, H. Kohno, and K. Miyake, J. Phys. Soc. Jpn. 
72 , 2914 (2003). 

[30] There Is a report that an s-wave superconductivity oc¬ 
curs In the Kondo lattice model within DMFT [32]. It 
employs the Infinite-dimensional Bethe lattice, and can¬ 
not be strictly compared with our results. 

[31] K. Masuda and D. Yamamoto, Phys. Rev. B 91, 104508 
(2015). 

[32] O. Bodensiek, R. Zitko, M. Vojta, M. Jarrell, and T. Pr- 
uschke, Phys. Rev. Lett. 110 , 146406 (2013). 


